Introducción a la Dendrocronología en R

Instrucciones

En cada sección podrás ejecutar y analizar los comandos básicos para leer, visualizar y realizar los análisis básicos utilizados en Dendrocronología. La mayoría de análisis que se llevaran acabo son iguales o muy similares a los que se realizan en los programas originales programados en FORTRAN llamados COFECHA y ARSTAN (https://www.ldeo.columbia.edu/tree-ring-laboratory/resources/software).

En este caso utilizaremos el ambiente de R y utilizaremos el paquete dplR programado y desarrollado por Andy Bunn y Mikko Korpela en el 2014 pueden bajar el articulo dándole clic al siguiente enlace. https://drive.google.com/uc?export=download&id=1ps7tZ0RTwuIr2P_h-ugHkbgAs3Jrdxaa

La idea y objetivo principal de estos ejercicios es que puedas familiarizarte con este tipo de ambiente, y puedas aprender haciendo. A lo largo de este documento hay un set de comandos y ejercicios que podrás ejecutar y modificar con tus propios archivos y o mediciones. Si lo modificaste y necesitas el archivo original lo puedes clonar de nuevo desde la pagina de github. https://github.com/paulszejner/Intro_dendrocrono.git

A continuación vamos a instalar y activar el dplR y se conectaran a la base de datos disponibles del international Tree Ring Data Bank (ITRDB) mediante la web

El primer paso de buenas costumbres y buenas practicas es iniciar con el ambiente de R limpio utilizando el comando rm(list=ls()) podrás limpiar todos los objetos que pueden estar activos en tu consola.

# para instalar el paquete de inertes hay que borrar el numeral que aparece antes del  "install.packages(dplR)" y ejecutar esa linea
# install.packages(dplR)

library(dplR) # librería para análisis dendrocronológicos

Warning: package ‘dplR’ was built under R version 3.6.2

Dat_ringwidth <-  read.rwl("https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/mexico/mexi042.rwl")

Attempting to automatically detect format.
Assuming a Tucson format file.
There appears to be a header in the rwl file
There are 20 series
1       RAN01A       1950    2006   0.001
2       RAN01B       1929    2006   0.001
3       RAN02A       1957    2006   0.001
4       RAN02B       1933    2006   0.001
5       RAN03A       1929    2006   0.001
6       RAN03B       1942    2006   0.001
7       RAN04A       1921    2006   0.001
8       RAN05A       1925    2006   0.001
9       RAN05B       1934    2006   0.001
10      RAN06A       1925    2006   0.001
11      RAN07A       1924    2006   0.001
12      RAN07B       1918    2006   0.001
13      RAN08A       1913    2006   0.001
14      RAN08B       1913    2006   0.001
15      RAN09A       1913    2006   0.001
16      RAN09B       1913    2006   0.001
17      RAN10A       1865    2006   0.001
18      RAN11B       1946    2006   0.001
19      RAN12A       1936    2006   0.001
20      RAN12B       1936    2006   0.001

Dat_Earlywood <- read.rwl("https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/mexico/mexi042e.rwl")

Attempting to automatically detect format.
Assuming a Tucson format file.
There appears to be a header in the rwl file
There are 20 series
1       RAN01A       1950    2006   0.001
2       RAN01B       1929    2006   0.001
3       RAN02A       1957    2006   0.001
4       RAN02B       1933    2006   0.001
5       RAN03A       1929    2006   0.001
6       RAN03B       1942    2006   0.001
7       RAN04A       1921    2006   0.001
8       RAN05A       1925    2006   0.001
9       RAN05B       1934    2006   0.001
10      RAN06A       1925    2006   0.001
11      RAN07A       1924    2006   0.001
12      RAN07B       1918    2006   0.001
13      RAN08A       1913    2006   0.001
14      RAN08B       1913    2006   0.001
15      RAN09A       1913    2006   0.001
16      RAN09B       1913    2006   0.001
17      RAN10A       1865    2006   0.001
18      RAN11B       1946    2006   0.001
19      RAN12A       1936    2006   0.001
20      RAN12B       1936    2006   0.001

Dat_latewood <- read.rwl("https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/mexico/mexi042l.rwl")

Attempting to automatically detect format.
Assuming a Tucson format file.
There appears to be a header in the rwl file
There are 17 series
1       RAN02A       1960    2006   0.001
2       RAN02B       1933    2006   0.001
3       RAN03A       1929    2006   0.001
4       RAN03B       1942    2006   0.001
5       RAN04A       1921    2006   0.001
6       RAN05A       1925    2006   0.001
7       RAN05B       1934    2006   0.001
8       RAN06A       1925    2006   0.001
9       RAN07A       1924    2006   0.001
10      RAN07B       1918    2006   0.001
11      RAN08A       1913    2006   0.001
12      RAN09A       1913    2006   0.001
13      RAN09B       1913    2006   0.001
14      RAN10A       1865    2006   0.001
15      RAN11B       1946    2006   0.001
16      RAN12A       1936    2006   0.001
17      RAN12B       1936    2006   0.001

# También se puede bajar  desde los FPTs del NOAA
# Dat_ftp <- read.rwl("ftp://ftp.ncdc.noaa.gov/pub/Data/paleo/treering/measurements/northamerica/mexico/mexi042.rwl")

str(Dat_ringwidth) # la función str te da la estructura del objeto llamado Datos

Classes ‘rwl’ and 'data.frame': 142 obs. of  20 variables:
 $ RAN01A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN01B: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN02A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN02B: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN03A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN03B: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN04A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN05A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN05B: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN06A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN07A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN07B: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN08A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN08B: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN09A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN09B: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN10A: num  1.08 1.47 1.7 1.08 1.48 ...
 $ RAN11B: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN12A: num  NA NA NA NA NA NA NA NA NA NA ...
 $ RAN12B: num  NA NA NA NA NA NA NA NA NA NA ...

### Grafica los datos para el control de calidad -------------------------------

plot.ts(Dat_ringwidth[,20])

plot.ts((Dat_ringwidth[,2]), col="darkred")

plot.ts((Dat_ringwidth[,3]), col="darkred")

plot.ts((Dat_ringwidth[,6]), col="darkred")

cuantas series de tiempo están en este archivo?

Podemos revisar todos los detalles del archivo con el cual vamos a trabajar

# cuantas series de tiempo están en este archivo? preguntando cuantas columnas  tenemos  es una forma de responder la pregunta
n_series <- ncol(Dat_ringwidth)
# ahora hacemos un plot con todas las serias  en su estado crudo o RAW  sin ningún tratamiento alguno
plot.ts(Dat_ringwidth[,20])
for(s in 1:n_series){lines(Dat_ringwidth[,s]) } # esta forma de  plotear todas las series es un for loop  usando cada columna  en cada uno de los loops

También podemos utilizar un comando especial del dplR llamado rwl.report()

rwl.report(Dat_ringwidth)

Number of dated series: 20 
Number of measurements: 1618 
Avg series length: 80.9 
Range:  142 
Span:  1865 - 2006 
Mean (Std dev) series intercorrelation: 0.4451521 (0.08958034)
Mean (Std dev) AR1: 0.7406 (0.08279931)
-------------
Years with absent rings listed by series 
    None 
-------------
Years with internal NA values listed by series 
    None 

En dplR hay otras opciones para graficar los datos que tenemos. Por ejemplo con mas orden y con algunas características que pueden ser útiles. La funcion que se utiliza es el famoso spaguetti plot spag.plot funcion del dplR

spag.plot(Dat_ringwidth, zfac = 1)

spag.plot(Dat_ringwidth, zfac = 3)

spag.plot(Dat_ringwidth, zfac = 0.3)

spag.plot(Dat_Earlywood, zfac = 0.3)

spag.plot(Dat_latewood, zfac = 0.3)

Como corroboramos el cofechado?

Luego de fechar visualmente las muestras de madera, se procede a su medición, Durante la medición se pueden cometer algunos errores, lo cual puede perjudicar la señal común de los árboles. Por lo tanto debemos corroborar las correlaciones entre todas las series antes de proceder a cualquier análisis.

#simple correlación entre series y series con su periodo en común
cor_matrix <- round(cor(Dat_ringwidth, use = "pairwise"), digits = 2) # la función cor  va  a utilizar una matriz o un data frame  para calcular  las correlaciones entre todas las columnas presentes en la base de datos. generando una matriz cuadrada  con le numero de columnas y  filas  igual al número de series.
cor_matrix

       RAN01A RAN01B RAN02A RAN02B RAN03A RAN03B RAN04A RAN05A RAN05B RAN06A RAN07A RAN07B
RAN01A   1.00   0.54   0.19   0.50   0.56   0.50   0.35   0.70   0.78   0.67   0.47   0.42
RAN01B   0.54   1.00   0.16   0.24   0.07   0.05   0.24   0.34   0.30   0.28  -0.18  -0.20
RAN02A   0.19   0.16   1.00   0.52   0.16   0.05   0.28   0.26   0.30   0.11   0.27   0.15
RAN02B   0.50   0.24   0.52   1.00   0.32   0.41   0.22   0.37   0.45   0.53   0.35   0.32
RAN03A   0.56   0.07   0.16   0.32   1.00   0.79   0.41   0.53   0.68   0.51   0.61   0.56
RAN03B   0.50   0.05   0.05   0.41   0.79   1.00   0.47   0.58   0.73   0.80   0.58   0.62
RAN04A   0.35   0.24   0.28   0.22   0.41   0.47   1.00   0.63   0.48   0.56   0.18   0.34
RAN05A   0.70   0.34   0.26   0.37   0.53   0.58   0.63   1.00   0.72   0.63   0.06   0.23
RAN05B   0.78   0.30   0.30   0.45   0.68   0.73   0.48   0.72   1.00   0.83   0.41   0.42
RAN06A   0.67   0.28   0.11   0.53   0.51   0.80   0.56   0.63   0.83   1.00   0.20   0.43
RAN07A   0.47  -0.18   0.27   0.35   0.61   0.58   0.18   0.06   0.41   0.20   1.00   0.66
RAN07B   0.42  -0.20   0.15   0.32   0.56   0.62   0.34   0.23   0.42   0.43   0.66   1.00
RAN08A   0.79   0.51   0.22   0.33   0.52   0.55   0.37   0.70   0.75   0.59  -0.04   0.25
RAN08B   0.80   0.51   0.28   0.40   0.58   0.61   0.51   0.76   0.82   0.68   0.05   0.28
RAN09A   0.60   0.09   0.17   0.55   0.55   0.66   0.58   0.50   0.60   0.67   0.45   0.62
RAN09B   0.56   0.07   0.15   0.44   0.62   0.79   0.61   0.59   0.63   0.69   0.38   0.65
RAN10A   0.74   0.35   0.24   0.35   0.58   0.61   0.34   0.71   0.79   0.60   0.10   0.31
RAN11B   0.33   0.24   0.44   0.31   0.01  -0.14  -0.32   0.13   0.14  -0.10  -0.02  -0.15
RAN12A   0.41   0.15   0.28   0.29   0.37   0.58   0.25   0.32   0.48   0.48   0.08   0.09
RAN12B   0.63   0.33   0.20   0.24   0.45   0.53   0.17   0.51   0.67   0.47   0.02   0.12
       RAN08A RAN08B RAN09A RAN09B RAN10A RAN11B RAN12A RAN12B
RAN01A   0.79   0.80   0.60   0.56   0.74   0.33   0.41   0.63
RAN01B   0.51   0.51   0.09   0.07   0.35   0.24   0.15   0.33
RAN02A   0.22   0.28   0.17   0.15   0.24   0.44   0.28   0.20
RAN02B   0.33   0.40   0.55   0.44   0.35   0.31   0.29   0.24
RAN03A   0.52   0.58   0.55   0.62   0.58   0.01   0.37   0.45
RAN03B   0.55   0.61   0.66   0.79   0.61  -0.14   0.58   0.53
RAN04A   0.37   0.51   0.58   0.61   0.34  -0.32   0.25   0.17
RAN05A   0.70   0.76   0.50   0.59   0.71   0.13   0.32   0.51
RAN05B   0.75   0.82   0.60   0.63   0.79   0.14   0.48   0.67
RAN06A   0.59   0.68   0.67   0.69   0.60  -0.10   0.48   0.47
RAN07A  -0.04   0.05   0.45   0.38   0.10  -0.02   0.08   0.02
RAN07B   0.25   0.28   0.62   0.65   0.31  -0.15   0.09   0.12
RAN08A   1.00   0.93   0.35   0.46   0.66   0.27   0.39   0.70
RAN08B   0.93   1.00   0.45   0.51   0.60   0.17   0.48   0.70
RAN09A   0.35   0.45   1.00   0.89   0.35  -0.27   0.24   0.17
RAN09B   0.46   0.51   0.89   1.00   0.51  -0.30   0.29   0.24
RAN10A   0.66   0.60   0.35   0.51   1.00   0.25   0.36   0.67
RAN11B   0.27   0.17  -0.27  -0.30   0.25   1.00  -0.04   0.30
RAN12A   0.39   0.48   0.24   0.29   0.36  -0.04   1.00   0.71
RAN12B   0.70   0.70   0.17   0.24   0.67   0.30   0.71   1.00

diag(cor_matrix)  #La función ¨diag¨ solamente extrae un el vector  diagonal  de  una matriz cuadrada,  y en este caso estamos sacando la correlación entre la misma serie. y por eso  deberíamos  obtener  el valor de 1 en toda la diagonal ya que estamos comparando lo mismo con lo mismo. 

RAN01A RAN01B RAN02A RAN02B RAN03A RAN03B RAN04A RAN05A RAN05B RAN06A RAN07A RAN07B RAN08A 
     1      1      1      1      1      1      1      1      1      1      1      1      1 
RAN08B RAN09A RAN09B RAN10A RAN11B RAN12A RAN12B 
     1      1      1      1      1      1      1 

Hay funciones en el dplR que pueden generar resultados muy similares a los de COFECHA, el cual es un programa que evalúa la sincronicidad de todas las series entre si. para hacer esta evaluación se dividen las series en segmentos correlativos y traslapados, para evaluar las correlaciones de todas las muestras por segmentos y sin tendencias que pueden sesgar los coeficientes de correlación.

A continuación se generaran unos ejemplos usando la función del dplR llamada corr.rwl.seg siempre es bueno poder utilizar los ejemplos que el mismo paquete provee. los cuales pueden acceder utilizando el siguiente comando, ?corr.rwl.seg en donde pueden leer todo sobre la función, y en la parte inferior de la pagina de ayuda verán algunos ejemplos que se pueden ejecutar en la consola.

En los siguientes ejemplos podrás ver que se ejecuta la misma función pero con diferentes argumentos. (los argumentos son los detalles que definirán como se evalúa y ejecuta la función). por ejemplo acá podemos aclarar el largo del segmento que se utiliza para comparar cada serie con el resto de las series seg.length=10, como también podemos indicarle que se remueva toda la autocorrelación que puede existir en las series prewhiten = T, generar un gráfico representando los resultados make.plot = T y delimitar el limite de confianza para su evaluación estadística pcrit = 0.1

nota: los colores rojos en las figuras indican los segmentos que presentan correlaciones por debajo del limite de confianza que se exije.

COFECHA_p_0.1 <- corr.rwl.seg(rwl = Dat_ringwidth, seg.length = 10, prewhiten = T, biweight = T, make.plot = T, pcrit = 0.1)

COFECHA_p_0.1 <- corr.rwl.seg(rwl = Dat_ringwidth, seg.length = 20, prewhiten = T, biweight = T, make.plot = T, pcrit = 0.1)

COFECHA_p_0.05 <- corr.rwl.seg(rwl = Dat_ringwidth, seg.length = 30, prewhiten = T, biweight = T, make.plot = T, pcrit = 0.05)

COFECHA_p_0.1 <- corr.rwl.seg(rwl = Dat_ringwidth, seg.length = 30, prewhiten = T, biweight = T, make.plot = T, pcrit = 0.1)

De todas las opciones que se pueden ejecutar, qué puedes concluir de este tipo de análisis? -que argumentos son importantes? -de que nos sirve este tipo de análisis?

COFECHA_p_0.1

$spearman.rho
       1900.1929 1915.1944 1930.1959 1945.1974 1960.1989  1975.2004
RAN01A        NA        NA        NA        NA 0.5052280 0.18620690
RAN01B        NA        NA 0.2556174 0.2142380 0.5924360 0.17775306
RAN02A        NA        NA        NA        NA 0.4896552 0.65250278
RAN02B        NA        NA        NA 0.4327030 0.4371524 0.38242492
RAN03A        NA        NA        NA 0.4967742 0.3147942 0.07897664
RAN03B        NA        NA        NA 0.5537264 0.4994438 0.36907675
RAN04A        NA        NA 0.4251390 0.4909900 0.6275862 0.66451613
RAN05A        NA        NA 0.1390434 0.2262514 0.4037820 0.39532814
RAN05B        NA        NA        NA 0.3619577 0.5688543 0.53058954
RAN06A        NA        NA 0.4233593 0.3971079 0.4518354 0.24983315
RAN07A        NA        NA 0.3294772 0.2840934 0.4371524 0.46963293
RAN07B        NA        NA 0.5715239 0.6907675 0.3810901 0.26318131
RAN08A        NA        NA 0.3023359 0.4745273 0.4665184 0.27919911
RAN08B        NA 0.5639600 0.3948832 0.3971079 0.6235818 0.52880979
RAN09A        NA        NA 0.5101224 0.4438265 0.3486096 0.38553949
RAN09B        NA        NA 0.6200222 0.6796440 0.5595106 0.57196885
RAN10A        NA 0.1350389 0.3005562 0.3241379 0.4527253 0.44115684
RAN11B        NA        NA        NA        NA 0.3788654 0.44382647
RAN12A        NA        NA        NA 0.3543938 0.3971079 0.49410456
RAN12B        NA        NA        NA 0.5648498 0.4936596 0.37931034

$p.val
       1900.1929    1915.1944    1930.1959    1945.1974    1960.1989    1975.2004
RAN01A        NA           NA           NA           NA 0.0024487314 1.615356e-01
RAN01B        NA           NA 0.0861259490 1.272391e-01 0.0003551341 1.729172e-01
RAN02A        NA           NA           NA           NA 0.0032976120 6.713042e-05
RAN02B        NA           NA           NA 8.863031e-03 0.0082473992 1.890727e-02
RAN03A        NA           NA           NA 2.882618e-03 0.0452618173 3.386263e-01
RAN03B        NA           NA           NA 8.888079e-04 0.0027389902 2.274730e-02
RAN04A        NA           NA 0.0099977278 3.216123e-03 0.0001395279 4.608887e-05
RAN05A        NA           NA 0.2310322792 1.141534e-01 0.0138720501 1.571371e-02
RAN05B        NA           NA           NA 2.503892e-02 0.0006288791 1.466418e-03
RAN06A        NA           NA 0.0102816158 1.531020e-02 0.0064640336 9.121239e-02
RAN07A        NA           NA 0.0379430229 6.402816e-02 0.0082473992 4.749131e-03
RAN07B        NA           NA 0.0005906781 1.913348e-05 0.0192657779 7.978581e-02
RAN08A        NA           NA 0.0522860139 4.351818e-03 0.0050176884 6.748935e-02
RAN08B        NA 0.0007045424 0.0158159438 1.531020e-02 0.0001560500 1.521954e-03
RAN09A        NA           NA 0.0022241898 7.391411e-03 0.0298326326 1.809184e-02
RAN09B        NA           NA 0.0001721577 2.804340e-05 0.0007800962 5.845124e-04
RAN10A        NA 0.2375856870 0.0533536507 4.048901e-02 0.0063673391 7.724380e-03
RAN11B        NA           NA           NA           NA 0.0198755491 7.391411e-03
RAN12A        NA           NA           NA 2.767294e-02 0.0153102020 3.032660e-03
RAN12B        NA           NA           NA 6.902227e-04 0.0030583065 1.975236e-02

$overall
             rho        p-val
RAN01A 0.4687628 1.571104e-04
RAN01B 0.2615542 1.091805e-02
RAN02A 0.4532840 7.664083e-04
RAN02B 0.4802881 1.240420e-05
RAN03A 0.4500615 2.722753e-05
RAN03B 0.5190476 7.485563e-06
RAN04A 0.5716716 1.240727e-08
RAN05A 0.3040876 2.980242e-03
RAN05B 0.4712274 2.077248e-05
RAN06A 0.4509711 1.461786e-05
RAN07A 0.3609079 5.009908e-04
RAN07B 0.4990158 6.170749e-07
RAN08A 0.3769783 2.426000e-04
RAN08B 0.5247908 6.316623e-08
RAN09A 0.4878802 6.511242e-07
RAN09B 0.6065775 0.000000e+00
RAN10A 0.3606800 2.242544e-04
RAN11B 0.3299825 5.507254e-03
RAN12A 0.4130697 2.124918e-04
RAN12B 0.5122037 3.975819e-06

$avg.seg.rho
1900.1929 1915.1944 1930.1959 1945.1974 1960.1989 1975.2004 
      NaN 0.3494994 0.3883709 0.4345351 0.4714794 0.3971969 

$flags
                RAN01A                 RAN01B                 RAN03A                 RAN05A 
           "1975.2004" "1945.1974, 1975.2004"            "1975.2004" "1930.1959, 1945.1974" 
                RAN10A 
           "1915.1944" 

$bins
     [,1] [,2]
[1,] 1900 1929
[2,] 1915 1944
[3,] 1930 1959
[4,] 1945 1974
[5,] 1960 1989
[6,] 1975 2004

$rwi
        RAN01A    RAN01B    RAN02A    RAN02B    RAN03A    RAN03B    RAN04A    RAN05A    RAN05B
1865        NA        NA        NA        NA        NA        NA        NA        NA        NA
1866        NA        NA        NA        NA        NA        NA        NA        NA        NA
1867        NA        NA        NA        NA        NA        NA        NA        NA        NA
1868        NA        NA        NA        NA        NA        NA        NA        NA        NA
1869        NA        NA        NA        NA        NA        NA        NA        NA        NA
1870        NA        NA        NA        NA        NA        NA        NA        NA        NA
1871        NA        NA        NA        NA        NA        NA        NA        NA        NA
1872        NA        NA        NA        NA        NA        NA        NA        NA        NA
1873        NA        NA        NA        NA        NA        NA        NA        NA        NA
1874        NA        NA        NA        NA        NA        NA        NA        NA        NA
1875        NA        NA        NA        NA        NA        NA        NA        NA        NA
1876        NA        NA        NA        NA        NA        NA        NA        NA        NA
1877        NA        NA        NA        NA        NA        NA        NA        NA        NA
1878        NA        NA        NA        NA        NA        NA        NA        NA        NA
1879        NA        NA        NA        NA        NA        NA        NA        NA        NA
1880        NA        NA        NA        NA        NA        NA        NA        NA        NA
1881        NA        NA        NA        NA        NA        NA        NA        NA        NA
1882        NA        NA        NA        NA        NA        NA        NA        NA        NA
1883        NA        NA        NA        NA        NA        NA        NA        NA        NA
1884        NA        NA        NA        NA        NA        NA        NA        NA        NA
1885        NA        NA        NA        NA        NA        NA        NA        NA        NA
1886        NA        NA        NA        NA        NA        NA        NA        NA        NA
1887        NA        NA        NA        NA        NA        NA        NA        NA        NA
1888        NA        NA        NA        NA        NA        NA        NA        NA        NA
1889        NA        NA        NA        NA        NA        NA        NA        NA        NA
1890        NA        NA        NA        NA        NA        NA        NA        NA        NA
1891        NA        NA        NA        NA        NA        NA        NA        NA        NA
1892        NA        NA        NA        NA        NA        NA        NA        NA        NA
1893        NA        NA        NA        NA        NA        NA        NA        NA        NA
1894        NA        NA        NA        NA        NA        NA        NA        NA        NA
1895        NA        NA        NA        NA        NA        NA        NA        NA        NA
1896        NA        NA        NA        NA        NA        NA        NA        NA        NA
1897        NA        NA        NA        NA        NA        NA        NA        NA        NA
1898        NA        NA        NA        NA        NA        NA        NA        NA        NA
1899        NA        NA        NA        NA        NA        NA        NA        NA        NA
1900        NA        NA        NA        NA        NA        NA        NA        NA        NA
1901        NA        NA        NA        NA        NA        NA        NA        NA        NA
1902        NA        NA        NA        NA        NA        NA        NA        NA        NA
1903        NA        NA        NA        NA        NA        NA        NA        NA        NA
1904        NA        NA        NA        NA        NA        NA        NA        NA        NA
1905        NA        NA        NA        NA        NA        NA        NA        NA        NA
1906        NA        NA        NA        NA        NA        NA        NA        NA        NA
1907        NA        NA        NA        NA        NA        NA        NA        NA        NA
1908        NA        NA        NA        NA        NA        NA        NA        NA        NA
1909        NA        NA        NA        NA        NA        NA        NA        NA        NA
1910        NA        NA        NA        NA        NA        NA        NA        NA        NA
1911        NA        NA        NA        NA        NA        NA        NA        NA        NA
1912        NA        NA        NA        NA        NA        NA        NA        NA        NA
1913        NA        NA        NA        NA        NA        NA        NA        NA        NA
1914        NA        NA        NA        NA        NA        NA        NA        NA        NA
         RAN06A    RAN07A      RAN07B    RAN08A    RAN08B     RAN09A      RAN09B      RAN10A
1865         NA        NA          NA        NA        NA         NA          NA          NA
1866         NA        NA          NA        NA        NA         NA          NA          NA
1867         NA        NA          NA        NA        NA         NA          NA          NA
1868         NA        NA          NA        NA        NA         NA          NA          NA
1869         NA        NA          NA        NA        NA         NA          NA          NA
1870         NA        NA          NA        NA        NA         NA          NA          NA
1871         NA        NA          NA        NA        NA         NA          NA          NA
1872         NA        NA          NA        NA        NA         NA          NA  0.94651776
1873         NA        NA          NA        NA        NA         NA          NA  0.90397343
1874         NA        NA          NA        NA        NA         NA          NA  0.57242195
1875         NA        NA          NA        NA        NA         NA          NA  1.12976698
1876         NA        NA          NA        NA        NA         NA          NA  0.76456242
1877         NA        NA          NA        NA        NA         NA          NA  1.14139152
1878         NA        NA          NA        NA        NA         NA          NA  0.72652422
1879         NA        NA          NA        NA        NA         NA          NA  1.05778896
1880         NA        NA          NA        NA        NA         NA          NA  1.01903887
1881         NA        NA          NA        NA        NA         NA          NA  1.33742461
1882         NA        NA          NA        NA        NA         NA          NA  0.75486804
1883         NA        NA          NA        NA        NA         NA          NA  1.10173068
1884         NA        NA          NA        NA        NA         NA          NA  0.83532752
1885         NA        NA          NA        NA        NA         NA          NA  0.84445707
1886         NA        NA          NA        NA        NA         NA          NA  1.07463069
1887         NA        NA          NA        NA        NA         NA          NA  0.80557260
1888         NA        NA          NA        NA        NA         NA          NA  1.08410699
1889         NA        NA          NA        NA        NA         NA          NA  0.85647449
1890         NA        NA          NA        NA        NA         NA          NA  1.14562403
1891         NA        NA          NA        NA        NA         NA          NA  0.60590956
1892         NA        NA          NA        NA        NA         NA          NA  1.13200724
1893         NA        NA          NA        NA        NA         NA          NA  0.85708224
1894         NA        NA          NA        NA        NA         NA          NA  0.84641965
1895         NA        NA          NA        NA        NA         NA          NA  1.19025571
1896         NA        NA          NA        NA        NA         NA          NA  0.96574328
1897         NA        NA          NA        NA        NA         NA          NA  1.41057017
1898         NA        NA          NA        NA        NA         NA          NA  1.07557579
1899         NA        NA          NA        NA        NA         NA          NA  0.76823276
1900         NA        NA          NA        NA        NA         NA          NA  1.27857764
1901         NA        NA          NA        NA        NA         NA          NA  0.75577425
1902         NA        NA          NA        NA        NA         NA          NA  1.31200133
1903         NA        NA          NA        NA        NA         NA          NA  0.71761940
1904         NA        NA          NA        NA        NA         NA          NA  0.85518087
1905         NA        NA          NA        NA        NA         NA          NA  0.88270042
1906         NA        NA          NA        NA        NA         NA          NA  0.92232665
1907         NA        NA          NA        NA        NA         NA          NA  0.98086993
1908         NA        NA          NA        NA        NA         NA          NA  1.19019025
1909         NA        NA          NA        NA        NA         NA          NA  0.83568290
1910         NA        NA          NA        NA        NA         NA          NA  0.89317761
1911         NA        NA          NA        NA        NA         NA          NA  1.13548355
1912         NA        NA          NA        NA        NA         NA          NA  0.89662268
1913         NA        NA          NA        NA        NA         NA          NA  0.95042806
1914         NA        NA          NA        NA        NA         NA          NA  0.92426289
        RAN11B    RAN12A    RAN12B
1865        NA        NA        NA
1866        NA        NA        NA
1867        NA        NA        NA
1868        NA        NA        NA
1869        NA        NA        NA
1870        NA        NA        NA
1871        NA        NA        NA
1872        NA        NA        NA
1873        NA        NA        NA
1874        NA        NA        NA
1875        NA        NA        NA
1876        NA        NA        NA
1877        NA        NA        NA
1878        NA        NA        NA
1879        NA        NA        NA
1880        NA        NA        NA
1881        NA        NA        NA
1882        NA        NA        NA
1883        NA        NA        NA
1884        NA        NA        NA
1885        NA        NA        NA
1886        NA        NA        NA
1887        NA        NA        NA
1888        NA        NA        NA
1889        NA        NA        NA
1890        NA        NA        NA
1891        NA        NA        NA
1892        NA        NA        NA
1893        NA        NA        NA
1894        NA        NA        NA
1895        NA        NA        NA
1896        NA        NA        NA
1897        NA        NA        NA
1898        NA        NA        NA
1899        NA        NA        NA
1900        NA        NA        NA
1901        NA        NA        NA
1902        NA        NA        NA
1903        NA        NA        NA
1904        NA        NA        NA
1905        NA        NA        NA
1906        NA        NA        NA
1907        NA        NA        NA
1908        NA        NA        NA
1909        NA        NA        NA
1910        NA        NA        NA
1911        NA        NA        NA
1912        NA        NA        NA
1913        NA        NA        NA
1914        NA        NA        NA
 [ reached getOption("max.print") -- omitted 92 rows ]

$seg.lag
[1] 15

$seg.length
[1] 30

$pcrit
[1] 0.1

$label.cex
[1] 1

attr(,"class")
[1] "list" "crs" 

Como hemos podido observar hay segmentos que no correlacionan con el conjunto de series. Por lo tanto podemos buscar los segmentos que no correlacionan con el resto de series utilizando el objeto resultado de la función corr.rwl.seg el cual se llama COFECHA_p_0.1y podemos indagar adentro de ese objeto utilizando COFECHA_p_0.1$flags.

COFECHA_p_0.1$flags

                RAN01A                 RAN01B                 RAN03A                 RAN05A 
           "1975.2004" "1945.1974, 1975.2004"            "1975.2004" "1930.1959, 1945.1974" 
                RAN10A 
           "1915.1944" 

#se puede especificar la muestra que te interesa corregir 
COFECHA_p_0.1$flags["RAN10A"] 

     RAN10A 
"1915.1944" 

COFECHA_p_0.1$flags["RAN05A"] 

                RAN05A 
"1930.1959, 1945.1974" 

# Graficar la muestra problemática con otras muestras que se ven bien

# En dplR los nombres de las columnas son las muestras y el nombre de las filas son los años
# Acá podemos extraer el año de inicio (el año mas antiguo del archivo  que tenemos)

Año_inicio <- as.numeric(rownames(Dat_ringwidth))[1]
time_series_convert <- ts(Dat_ringwidth[,c("RAN08B","RAN10A")],start = Año_inicio, frequency =1 )
plot(time_series_convert)

# Selección de un solo segmento de las series
segmento_Problemas <- window(time_series_convert, start=1915, end=1944) 
segmento_SIN_Problemas <- window(time_series_convert,1970,1990) 

# visualizamos el segmento con problemas con respecto a una seria sin problemas
plot(segmento_Problemas[,"RAN08B"])
lines(segmento_Problemas[,"RAN10A"], lwd=2)

# Tambien se puede identificar la columna por su numero
plot(segmento_SIN_Problemas[,1]) #"RAN08B"
lines(segmento_SIN_Problemas[,2], lwd=2) #"RAN10A"

# Ajustamos los  limites del Eje "Y"
plot(segmento_SIN_Problemas[,1], ylim= range(segmento_SIN_Problemas))
lines(segmento_SIN_Problemas[,2], lwd=2)

Remoción de tendencias (detrending), estandarizacion y extraccion de la señal buscada

Ejemplo ficticio

Primero podremos fabricar una serie de tiempo que se apegue a la teoría de el modelo linear agregado “Linear aggregated model” en donde se suman distintas fuentes de variabilidad como el error estocástico, random, Disturbios endógenos y exógenos, Variabilidad climática tendencias relacionadas con la edad.

Si quieres una buena presentación de lo que es el modelo agregado de crecimiento de los anillos de los árboles has clic!

https://drive.google.com/uc?export=download&id=1yk4RO9h6o9Xd-dbryFErCsQ302xW-9mf

# Acá vamos a generar y crear  diferentes vectores de  de 200 años, con diferentes  magnitudes, tendencias  y variabilidad. 

Edad <- 1*(1/1.04)^c(1:200) # formula de una curva exponencial negativa
clima <- rnorm(200,mean= 1, sd = 0.3) # se generara un vector con variabilidad  aleatoria con media de uno y desviación de 0.3
disturbios <- c(rep(x = 0, 100), seq(from = 0.8,to = 0.001, length.out = 30),rep(0,70))  #disturbio agregado en el año 101 con efecto de 30 años
Error <- rnorm(200, mean = 0.3,sd = 0.1 )

# luego sumamos cada uno de estos factores. 

  Agregado <- Error+disturbios+clima+Edad

El layout es útil para hacer plots en la misma ventana esta basado en una matriz en donde los números correlativos 1…x es el orden y posición en que hacemos el plot

layout(matrix(1:2, nrow = 2, ncol = 1),widths = 3, heights = c(2,2), respect = T )
par(mar=c(2,4,3,2))
plot.ts(Error,  ylim=c(0,2), ylab="variables")
lines(disturbios)

lines(clima)
lines(Edad)

plot.ts(Agregado, main="variabilidad agregada 
Error+clima+edad+disturbio")

NA
NA

#Otro ejmplo de uso del layout es:

layout(matrix(c(1,0,
                2,5,
                3,5,
                4,0), nrow = 4, ncol = 2, byrow = T),widths = c(2,3), heights = c(1,1), respect = T )

layout.show(5) # esto solo  enseña cual es la configuración de los 5 plots.

par(mar=c(1,4,1,1))

plot.ts(Error,  ylim=c(0,2))
plot.ts(disturbios, ylim=c(0,2))

plot.ts(clima, col="darkblue", ylim=c(0,2))
plot.ts(Edad,col="darkgreen", lwd=3, ylim=c(0,2))

par(mar=c(2,4,3,2))
plot.ts(Agregado, main="variabilidad agregada
Error+clima+edad+disturbio")

ahora vemos el effecto de la extraccion de la señal climatica y la edad utilizando el detrending o remocion de tendencias

Primero revisamos la correlación que excite en nuestra seria agregada versusla variabilidad del clima

plot(Agregado,clima) # señal del "agregado" versus el clima

cor(Agregado,clima) # señal del "agregado" versus el clima

[1] 0.7850268

Como podemos observar esta relación contiene variabilidad ajena al clima. Recordemos que el agregado contiene error random, un disturbio, y el efecto de la edad, por lo que la correlación no es perfecta y se puede mejorar utilizando algunas tecnicas que a continuación se realizaran.

En los siguientes pasos vamos a experimentar algunas técnicas para remover el efecto de la edad, y el efecto del disturbio utilizando herramientas de estadística y de moldeamiento linear.

Las funciones que se pueden utilizar para esta remoción de tendencias son: detrend.series la que funciona para remover la tendencia a series en especifico y detrend para remover la tendencia simultáneamente a muchas series.

detrend.series(Agregado, method = "ModNegExp") # en este caso  utilizamos  la función  con el argumento del método el cual  se utiliza para el ajuste linear. y se selecciona  el Modo exponencial negativo para que el programa  ajuste una función exponencial negativa  a la serie del agregado y se puedan estimar los residuales de este ajuste, para luego estandarizar con una media de uno.

  [1] 0.8892685 1.0589656 0.9360061 0.8989878 1.0806626 0.7616235 0.8659676 1.1813124 1.0807743
 [10] 1.3667037 0.9684134 0.7472339 0.8841024 0.9489908 1.1107495 0.8847535 1.3900642 1.2015906
 [19] 1.2880488 0.7665305 1.5358611 1.3718816 0.7380057 1.1608772 0.9244117 0.9109595 1.0217263
 [28] 0.8784879 0.7531410 1.0960533 0.7258404 1.2166896 1.0313231 1.1310552 1.0618982 1.1537623
 [37] 0.4371164 0.9975767 0.9425944 1.0080886 0.6383375 1.2119826 1.1637763 1.2937458 0.5726399
 [46] 0.6182085 0.6953434 0.7896170 0.5554152 1.1392820 0.8999486 0.9766777 0.9436739 0.6657177
 [55] 1.0106053 1.0762284 0.9536822 1.2657869 1.1920575 0.6836828 1.3282376 1.2567639 1.0418295
 [64] 1.0001245 0.8078312 0.3458015 1.1636359 0.9013626 0.9512976 0.7632315 0.6631120 1.1954283
 [73] 1.1767780 1.4005029 1.1481960 0.8678046 1.1834661 0.6788073 0.7486638 0.9639214 0.8654230
 [82] 1.1563842 0.6309422 0.7019131 0.8455205 1.0019664 1.0576932 0.9383308 1.1203158 1.2628665
 [91] 0.6804501 0.7941744 0.7239737 0.8720152 1.1986117 1.0364070 1.2734452 0.9738866 1.0383842
[100] 1.0628808 1.5075546 1.5111480 1.1580627 1.5885639 1.1856912 1.5396366 1.1937287 1.3932143
[109] 1.6262490 1.1891143 1.0530765 1.6135760 1.4230432 0.9393754 1.3587604 1.2656352 1.8178658
[118] 1.2019795 1.2200660 0.8985697 0.9087693 0.9660080 1.5307006 1.6025178 1.2778207 1.0682771
[127] 1.4654616 1.1397496 0.7996813 1.1061921 0.8661061 1.1777366 1.0745332 1.0871287 0.7761246
[136] 1.6282523 1.0232410 0.7490499 0.8115941 1.0399866 1.0614687 0.8280863 0.8581606 0.8293600
[145] 0.6457492 0.9387977 0.4966882 0.7864802 0.7471034 0.7880427 0.4957430 0.5279081 0.9618160
[154] 0.9011958 0.6633285 1.0715933 1.2731971 1.0533998 1.1854227 1.0138245 0.6997306 0.7084341
[163] 1.1790340 0.3957535 0.8784445 0.8265876 0.8351693 0.9660597 0.4892477 1.1636636 1.0742514
[172] 0.8718536 0.7510640 1.3524006 1.1944380 1.1392958 1.4362231 1.0227383 1.1815778 1.1293985
[181] 0.5789783 0.9879183 0.7997471 0.9424154 0.7115777 0.9711359 0.6414401 1.0677674 0.7226057
[190] 1.3494602 0.7957381 0.7176207 1.0862772 1.0312411 0.8419485 0.6633071 0.9329478 0.8196898
[199] 0.4149680 1.1327574

detrend.series(Agregado, method = "Spline", nyrs = 20) # en ese caso se ajusta una curva mas flexible  llamada Qubic spline  con 40 años de flexibilidad.

  [1] 0.9528747 1.1230354 0.9831914 0.9355689 1.1139465 0.7772158 0.8739135 1.1787941 1.0682022
 [10] 1.3412858 0.9456065 0.7254105 0.8504683 0.9011375 1.0383936 0.8137412 1.2596633 1.0774650
 [19] 1.1493790 0.6847256 1.3812206 1.2505005 0.6856336 1.1018792 0.8968849 0.9019225 1.0288310
 [28] 0.8957424 0.7735307 1.1281845 0.7462463 1.2474646 1.0558624 1.1596281 1.0941425 1.1981953
 [37] 0.4580241 1.0531129 1.0015537 1.0782826 0.6880754 1.3192289 1.2849438 1.4563211 0.6588235
 [46] 0.7245095 0.8236984 0.9364868 0.6533244 1.3188261 1.0208796 1.0831567 1.0218342 0.7031859
 [55] 1.0409582 1.0830845 0.9415337 1.2328322 1.1534516 0.6617484 1.2937654 1.2402286 1.0477336
 [64] 1.0278111 0.8472114 0.3679272 1.2453498 0.9633636 1.0091716 0.7995354 0.6833905 1.2112058
 [73] 1.1771636 1.3939720 1.1477418 0.8780891 1.2184826 0.7130079 0.8015111 1.0487667 0.9538378
 [82] 1.2866232 0.7058986 0.7852132 0.9405496 1.1045934 1.1548018 1.0157414 1.2043851 1.3503228
 [91] 0.7232773 0.8352957 0.7481497 0.8788676 1.1717098 0.9798715 1.1628948 0.8586945 0.8838114
[100] 0.8738817 1.2004067 1.1710326 0.8783236 1.1853974 0.8749825 1.1284736 0.8723944 1.0183528
[109] 1.1923437 0.8767336 0.7816166 1.2061272 1.0723862 0.7141116 1.0417716 0.9793101 1.4220728
[118] 0.9528195 0.9802513 0.7299744 0.7431628 0.7918251 1.2556314 1.3187715 1.0599640 0.8972543
[127] 1.2507093 0.9911330 0.7088452 0.9974694 0.7923820 1.0906042 1.0062532 1.0299619 0.7450785
[136] 1.5885953 1.0198647 0.7657834 0.8528242 1.1251329 1.1852628 0.9567225 1.0268887 1.0271326
[145] 0.8253311 1.2318386 0.6644463 1.0626131 1.0093472 1.0538240 0.6494902 0.6709563 1.1780991
[154] 1.0623496 0.7544135 1.1823669 1.3758035 1.1280193 1.2724910 1.1015809 0.7741625 0.7990919
[163] 1.3530791 0.4600255 1.0259653 0.9616069 0.9597523 1.0884437 0.5369976 1.2380649 1.1064015
[172] 0.8700865 0.7278903 1.2783946 1.1099035 1.0507383 1.3283601 0.9585141 1.1316721 1.1122647
[181] 0.5878879 1.0326856 0.8575007 1.0310883 0.7896792 1.0862922 0.7193635 1.1957944 0.8070821
[190] 1.5046210 0.8886120 0.8049052 1.2269967 1.1772235 0.9743401 0.7789722 1.1109188 0.9881918
[199] 0.5052127 1.3887114

# Acá repetimos lo mismo  pero le asignamos un nombre a cada  objeto nuevo
det_exp <- detrend.series(Agregado, method = "ModNegExp")

det_sp <- detrend.series(Agregado, method = "Spline", nyrs = 20)

layout(matrix(1:2, nrow = 2, ncol = 1),widths = 3, heights = c(2,2), respect = T )

par(mar=c(2,4,3,2))

plot.ts(det_exp)
lines(det_sp, col="red")

abline(h=1)

plot(det_exp,clima)

points(det_sp,clima, col="red", pch=19)
text(x = 0.6, y = 1.5,labels = round(cor(det_exp,clima),2))

text(x = 0.6, y = 1,labels = round(cor(det_sp,clima),2), col= "red")

cor(det_sp,Error)

[1] 0.3455075

Conclusión.

El spline de 25 años logra minimizar el efecto de la edad y de ese único evento de disturbio.

Otro ejemplo del modelo agregado es cuando observamos como el CO2 atmosférico o el dióxido de carbono. el CO2 varia mensualmente a lo largo de la historia. ver https://gml.noaa.gov/ccgg/trends/

require(graphics)
m <- decompose(co2)
plot(m)

Desarrollo de la cronología maestra

Luego de remover las tendencias ajenas al clima procedemos a generar la cronología maestra utilizando todas las series correctamente fechadas. para generar la cronología media de anchos de anillos podemos utilizar la función llamada chron Esta función construye una cronología del valor medio de todas las serias previamente estandarizadas. Típicamente a partir de un data.frame de anchos de anillo sin tendencia según lo producido por detrend.

# Recordamos  que tipo de  datos  tenemos  
rwl.report(Dat_ringwidth)

Number of dated series: 20 
Number of measurements: 1618 
Avg series length: 80.9 
Range:  142 
Span:  1865 - 2006 
Mean (Std dev) series intercorrelation: 0.4451521 (0.08958034)
Mean (Std dev) AR1: 0.7406 (0.08279931)
-------------
Years with absent rings listed by series 
    None 
-------------
Years with internal NA values listed by series 
    None 

# Podemos explotar las estadísticas descriptivas  de la base de datos
rwl.stats(Dat_ringwidth)

# Luego pasamos a remover las tendencias de la edad y de posibles disturbios 
Dat_ringwidth_index <- detrend(Dat_ringwidth,method = "Spline", nyrs = 50, make.plot = T)

NA

Después de remover las tendencias y estandarizar las series medidas, podemos calcular sus estadísticos para evaluar la señal común a lo largo de la cronologista.

Statistics <- rwi.stats.running (Dat_ringwidth_index, window.length = 20)

Dat_ringwidth_cronologia <- chron(Dat_ringwidth_index, biweight = T,prewhiten = T)

plot(Dat_ringwidth_cronologia)